1. Set up

#.........................Load Libraries.........................
library(janitor)
library(tidyverse)
library(ggeffects)
library(sjlabelled)
library(lme4)
library(DHARMa)
library(parameters)
library(effects)
library(MuMIn)
library(waffle)
library(showtext)
library(lmerTest)
library(here)

#.........................Load Fonts.........................
font_add_google(name = "Marcellus", family = "marc")
font_add_google(name = "Josefin Sans", family = "josefin")
font_add_google(name = "Lexend", family = "lex")

#...............Assigning a Standard ggplot Theme............

# ggplot theme
theme_set(theme_bw() +
            theme(panel.grid = element_blank(), 
                  plot.background = element_blank()))

# abline color
ref_line_col <- "grey"

# model prediction color
model_col <- "darkblue"

#.......Reading in Data for all Seven Traits and Treatments.........
 
# Each data frame is read in, then specimen_id column is separated to create new 
# columns, then we create an individual ID column from date, site, and individual.

# using read_csv and `here` package to get rid of file paths
height_all <- read_csv(here("data", "preliminary_data", "final_height.csv")) |> 
  # separate specimen_id column by - to create new date, site, individual, and specimen columns
  separate_wider_delim(specimen_id, 
                       delim = "-", 
                       names = c("date", "site", "individual", "specimen"), 
                       cols_remove = FALSE) |> 
  # create an individual id column
  unite("individual_id", date:individual, remove = TRUE, sep = "-") 


width_all <- read_csv(here("data", "preliminary_data", "final_width.csv")) |> 
  separate_wider_delim(specimen_id,
                       delim = "-",                                                                     names = c("date", "site", "individual", "specimen"),
                       cols_remove = FALSE) |> 
  unite("individual_id", date:individual, remove = TRUE, sep = "-")

perimeter_all <- read_csv(here("data", "preliminary_data", "final_perimeter.csv")) |> 
  separate_wider_delim(specimen_id, 
                       delim = "-", 
                       names = c("date", "site", "individual", "specimen"),
                       cols_remove = FALSE) |> 
  unite("individual_id", date:individual, remove = TRUE, sep = "-")

surf_all <- read_csv(here("data", "preliminary_data", "final_surf.csv")) |> 
  separate_wider_delim(specimen_id, 
                       delim = "-", 
                       names = c("date", "site", "individual", "specimen"), 
                       cols_remove = FALSE) |> 
  unite("individual_id", date:individual, remove = TRUE, sep = "-")

tdmc_all <- read_csv(here("data", "preliminary_data", "final_tdmc.csv")) |> 
  separate_wider_delim(specimen_id, 
                       delim = "-", 
                       names = c("date", "site", "individual", "specimen"), 
                       cols_remove = FALSE) |> 
  unite("individual_id", date:individual, remove = TRUE, sep = "-")

thickness_all <- read_csv(here("data", "preliminary_data", "final_thickness.csv")) |> 
  separate_wider_delim(specimen_id, 
                       delim = "-", 
                       names = c("date", "site", "individual", "specimen"),
                       cols_remove = FALSE) |> 
  unite("individual_id", date:individual, remove = TRUE, sep = "-")

volume_all <- read_csv(here("data", "preliminary_data", "final_volume.csv")) |>
  separate_wider_delim(specimen_id, 
                       delim = "-", 
                       names = c("date", "site", "individual", "specimen"), 
                       cols_remove = FALSE) |> 
  unite("individual_id", date:individual, remove = TRUE, sep = "-")

2. Linear mixed effect models

a. initial values predicted by rehydrate values and/or species

This section includes code to fit linear mixed effect models to answer the question: how well do rehydrated trait values predict initial trait values? Each trait is predicted by two models: initial values as a function of rehydrated values with a species interaction, or intiial values as a function of rehydrated values without a species interaction - both models have individual_id included as a random effect. We view model summaries using summary().

Height

Width

#..........................Width Models.........................

width_reg_ir <- lmer(width_i ~ width_r * species + (1 | individual_id), data = width_all)
summary(width_reg_ir)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_r * species + (1 | individual_id)
##    Data: width_all
## 
## REML criterion at convergence: 16.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8108 -0.3944 -0.1178  0.3648  5.3792 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.00000  0.0000  
##  Residual                  0.04315  0.2077  
## Number of obs: 162, groups:  individual_id, 54
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           0.00229    0.11672 140.00000   0.020   0.9844    
## width_r               1.03003    0.03538 140.00000  29.115   <2e-16 ***
## speciesCC             0.02795    0.15960 140.00000   0.175   0.8613    
## speciesCO             0.03873    0.22403 140.00000   0.173   0.8630    
## speciesCYOS           0.16350    0.16650 140.00000   0.982   0.3278    
## speciesDL             0.17353    0.77701 140.00000   0.223   0.8236    
## speciesDP             0.03608    0.13583 140.00000   0.266   0.7909    
## speciesEGME           0.47356    0.19897 140.00000   2.380   0.0187 *  
## speciesGR            -0.12681    0.25299 140.00000  -0.501   0.6170    
## speciesNAND           0.07054    0.16759 140.00000   0.421   0.6745    
## speciesPOLA           0.12897    0.67780 140.00000   0.190   0.8494    
## speciesR              0.04267    0.14411 140.00000   0.296   0.7676    
## width_r:speciesCC    -0.03103    0.03677 140.00000  -0.844   0.4002    
## width_r:speciesCO    -0.02300    0.06060 140.00000  -0.379   0.7049    
## width_r:speciesCYOS  -0.03545    0.04123 140.00000  -0.860   0.3914    
## width_r:speciesDL    -0.18023    0.65790 140.00000  -0.274   0.7845    
## width_r:speciesDP    -0.01653    0.04568 140.00000  -0.362   0.7180    
## width_r:speciesEGME  -0.15968    0.06262 140.00000  -2.550   0.0118 *  
## width_r:speciesGR     0.07176    0.10347 140.00000   0.694   0.4891    
## width_r:speciesNAND  -0.03328    0.04590 140.00000  -0.725   0.4695    
## width_r:speciesPOLA  -0.06468    0.19565 140.00000  -0.331   0.7415    
## width_r:speciesR     -0.05309    0.03821 140.00000  -1.389   0.1669    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
width_reg_ir_nospecies <- lmer(width_i ~ width_r + (1 | individual_id), data = width_all)
summary(width_reg_ir_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_r + (1 | individual_id)
##    Data: width_all
## 
## REML criterion at convergence: -38.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0894 -0.3833 -0.1711  0.4370  5.3312 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.00000  0.0000  
##  Residual                  0.04256  0.2063  
## Number of obs: 162, groups:  individual_id, 54
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 8.645e-02  2.502e-02 1.600e+02   3.456 0.000703 ***
## width_r     9.939e-01  4.949e-03 1.600e+02 200.819  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## width_r -0.762
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')

### Perimeter

perimeter_reg_ir <- lmer(per_i ~ per_r * species + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ir)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_r * species + (1 | individual_id)
##    Data: perimeter_all
## 
## REML criterion at convergence: 1298.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4623 -0.1754 -0.0524  0.1051  5.7167 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept)  95.52    9.774  
##  Residual                  317.46   17.817  
## Number of obs: 153, groups:  individual_id, 52
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)        11.94377    9.72785  88.63228   1.228   0.2228    
## per_r               0.89901    0.06292 128.62071  14.289  < 2e-16 ***
## speciesCC         -12.00315   14.14248 102.00347  -0.849   0.3980    
## speciesCO         -13.67406   16.07381 104.35426  -0.851   0.3969    
## speciesCYOS        -6.27899   15.71169 113.52599  -0.400   0.6902    
## speciesDP          -1.15338   13.35059  94.16302  -0.086   0.9313    
## speciesEGME       -10.09572   14.24611 108.75008  -0.709   0.4800    
## speciesGR          -7.51742   16.61084  85.79641  -0.453   0.6520    
## speciesNAND       -29.82349   34.36413  87.86417  -0.868   0.3878    
## speciesPOLA       188.18674   27.83589 131.21107   6.761 4.09e-10 ***
## speciesR          -10.42931   13.45056  89.53546  -0.775   0.4402    
## per_r:speciesCC     0.12080    0.13574 129.10913   0.890   0.3752    
## per_r:speciesCO     0.13346    0.13865 132.86649   0.963   0.3375    
## per_r:speciesCYOS   0.10090    0.08751 132.78026   1.153   0.2510    
## per_r:speciesDP     0.17766    0.07971 126.37640   2.229   0.0276 *  
## per_r:speciesEGME   0.06577    0.20012 107.37549   0.329   0.7431    
## per_r:speciesGR     0.09556    0.11385 125.76078   0.839   0.4029    
## per_r:speciesNAND   0.35021    0.22788 107.56440   1.537   0.1273    
## per_r:speciesPOLA  -1.04870    0.13015 110.28003  -8.058 1.01e-12 ***
## per_r:speciesR      0.09930    0.10495 131.89298   0.946   0.3458    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perimeter_reg_ir_nospecies <- lmer(per_i ~ per_r + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ir_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_r + (1 | individual_id)
##    Data: perimeter_all
## 
## REML criterion at convergence: 1434.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.3939 -0.2103 -0.1346  0.0058  4.3400 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept)  72.39    8.508  
##  Residual                  627.45   25.049  
## Number of obs: 153, groups:  individual_id, 52
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   5.56937    4.08128  67.48997   1.365    0.177    
## per_r         1.00652    0.02865 104.53362  35.138   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## per_r -0.814

Surface Area

surf_reg_ir <- lmer(surf_i ~ surf_r * species + (1 | individual_id), data = surf_all)
summary(surf_reg_ir)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_r * species + (1 | individual_id)
##    Data: surf_all
## 
## REML criterion at convergence: 1188.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1047 -0.1974  0.0000  0.1264  7.5343 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept)  81.54    9.03   
##  Residual                  117.26   10.83   
## Number of obs: 155, groups:  individual_id, 52
## 
## Fixed effects:
##                     Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         -1.58126    7.81026 107.08809  -0.202   0.8399    
## surf_r               1.12908    0.08684 134.94728  13.002   <2e-16 ***
## speciesCC           15.42369    9.77328  92.22910   1.578   0.1180    
## speciesCO            0.89677   10.38244 105.62177   0.086   0.9313    
## speciesCYOS          3.29739   11.49479 122.93006   0.287   0.7747    
## speciesDP            9.75021    9.20383 103.96275   1.059   0.2919    
## speciesEGME          2.18597    9.97535  97.72563   0.219   0.8270    
## speciesGR           -1.55355   15.32926  82.71308  -0.101   0.9195    
## speciesNAND         -3.82483   16.69643 128.64243  -0.229   0.8192    
## speciesPOLA         -4.48724   20.82619 132.22363  -0.215   0.8297    
## speciesR             0.15977   10.70600  82.79907   0.015   0.9881    
## surf_r:speciesCC    -0.16124    0.08978 134.80058  -1.796   0.0747 .  
## surf_r:speciesCO     0.05879    0.48435 128.38318   0.121   0.9036    
## surf_r:speciesCYOS   0.01105    0.13090 120.77518   0.084   0.9329    
## surf_r:speciesDP     0.10850    0.09911 134.94484   1.095   0.2755    
## surf_r:speciesEGME  -0.10908    0.12497 123.28075  -0.873   0.3844    
## surf_r:speciesGR     0.01837    0.17066 111.20076   0.108   0.9145    
## surf_r:speciesNAND   0.25954    0.34672 117.31254   0.749   0.4556    
## surf_r:speciesPOLA  -0.04449    0.17205 109.89074  -0.259   0.7964    
## surf_r:speciesR     -0.02750    0.15965  94.00692  -0.172   0.8636    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
surf_reg_ir_nospecies <- lmer(surf_i ~ surf_r + (1 | individual_id), data = surf_all)
summary(surf_reg_ir_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_r + (1 | individual_id)
##    Data: surf_all
## 
## REML criterion at convergence: 1273
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2860 -0.2405 -0.0934  0.1431  7.0725 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 164.2    12.82   
##  Residual                  128.3    11.33   
## Number of obs: 155, groups:  individual_id, 52
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   8.15738    2.38648  80.79972   3.418  0.00099 ***
## surf_r        1.03768    0.01904 145.54351  54.505  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## surf_r -0.534

Weight

tdmc_reg_ir <- lmer(weight_i ~ weight_r * species + (1 | individual_id), data = tdmc_all)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(tdmc_reg_ir)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_r * species + (1 | individual_id)
##    Data: tdmc_all
## 
## REML criterion at convergence: 2556.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2441 -0.2221 -0.0195  0.2035  7.5446 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept)       0     0    
##  Residual                  1382191  1176    
## Number of obs: 159, groups:  individual_id, 55
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)          -859.71069  389.57350  139.00000  -2.207  0.02897 *  
## weight_r                1.80897    0.03927  139.00000  46.061  < 2e-16 ***
## speciesCC            1948.48497  582.33898  139.00000   3.346  0.00105 ** 
## speciesCO             796.50963  592.86217  139.00000   1.343  0.18130    
## speciesCYOS          1548.97927  816.97129  139.00000   1.896  0.06004 .  
## speciesDP             826.85381  515.37765  139.00000   1.604  0.11090    
## speciesEGME           983.09662  599.62798  139.00000   1.640  0.10337    
## speciesGR            -290.19721  917.05077  139.00000  -0.316  0.75214    
## speciesNAND           898.08911 1658.39242  139.00000   0.542  0.58900    
## speciesPOLA           672.30438 1364.35381  139.00000   0.493  0.62296    
## speciesR              885.35159  604.43919  139.00000   1.465  0.14525    
## weight_r:speciesCC     -0.81629    0.05323  139.00000 -15.335  < 2e-16 ***
## weight_r:speciesCO     -0.50096    0.55557  139.00000  -0.902  0.36878    
## weight_r:speciesCYOS    0.61326    0.34120  139.00000   1.797  0.07445 .  
## weight_r:speciesDP      0.74481    0.17310  139.00000   4.303 3.16e-05 ***
## weight_r:speciesEGME   -0.30748    0.58338  139.00000  -0.527  0.59898    
## weight_r:speciesGR     -0.08245    0.09342  139.00000  -0.883  0.37901    
## weight_r:speciesNAND   -0.27214    1.17642  139.00000  -0.231  0.81740    
## weight_r:speciesPOLA   -0.35454    0.62942  139.00000  -0.563  0.57415    
## weight_r:speciesR      -0.44324    0.23046  139.00000  -1.923  0.05649 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
tdmc_reg_ir_nospecies <- lmer(weight_i ~ weight_r + (1 | individual_id), data = tdmc_all)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(tdmc_reg_ir_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_r + (1 | individual_id)
##    Data: tdmc_all
## 
## REML criterion at convergence: 2907
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0646 -0.1959 -0.1068  0.1539  6.4238 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 2316648  1522    
##  Residual                  3891255  1973    
## Number of obs: 159, groups:  individual_id, 55
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept) 832.53896  295.38612  69.63116   2.818  0.00628 ** 
## weight_r      1.39800    0.04207 153.61803  33.228  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## weight_r -0.459
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

Thickness

thick_reg_ir <- lmer(thick_i ~ thick_r * species + (1 | individual_id), data = thickness_all)
summary(thick_reg_ir)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_r * species + (1 | individual_id)
##    Data: thickness_all
## 
## REML criterion at convergence: -96.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2260 -0.1919 -0.0432  0.0915  7.2546 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.04666  0.2160  
##  Residual                  0.01759  0.1326  
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)           0.120542   0.114638  73.651818   1.052   0.2965  
## thick_r               0.139866   0.320579 102.193603   0.436   0.6635  
## speciesCC             0.010940   0.220029 122.641583   0.050   0.9604  
## speciesCO             0.363711   0.331424 132.362361   1.097   0.2745  
## speciesCYOS           0.457297   0.237785 136.671480   1.923   0.0565 .
## speciesDL             0.049458   0.459266 139.979514   0.108   0.9144  
## speciesDP             0.124849   0.222560 127.567063   0.561   0.5758  
## speciesEGME           0.001116   0.220785 123.673763   0.005   0.9960  
## speciesGR             0.360171   1.022235 108.311223   0.352   0.7253  
## speciesNAND           0.059782   0.373055 139.094404   0.160   0.8729  
## speciesPOLA           0.045172   1.568954  99.927148   0.029   0.9771  
## speciesR             -0.102688   0.194661  68.868970  -0.528   0.5995  
## thick_r:speciesCC     0.616505   0.480729 108.357418   1.282   0.2024  
## thick_r:speciesCO    -0.029348   0.664557 101.790886  -0.044   0.9649  
## thick_r:speciesCYOS   0.198255   0.886893  97.868959   0.224   0.8236  
## thick_r:speciesDL    -0.139866   2.330202  94.656268  -0.060   0.9523  
## thick_r:speciesDP     0.495902   1.270180 135.929443   0.390   0.6968  
## thick_r:speciesEGME   0.281383   1.275901 110.142813   0.221   0.8259  
## thick_r:speciesGR     0.673354   0.797153 126.255219   0.845   0.3999  
## thick_r:speciesNAND  -0.263202   1.800285  94.824969  -0.146   0.8841  
## thick_r:speciesPOLA  -0.497009  14.474206  94.516655  -0.034   0.9727  
## thick_r:speciesR      0.765868   0.655984 126.848214   1.168   0.2452  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
thick_reg_ir_nospecies <- lmer(thick_i ~ thick_r + (1 | individual_id), data = thickness_all)
summary(thick_reg_ir_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_r + (1 | individual_id)
##    Data: thickness_all
## 
## REML criterion at convergence: -63.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.2443 -0.2138 -0.0588  0.1191  7.2879 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.05910  0.2431  
##  Residual                  0.01742  0.1320  
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.13422    0.04655 56.48883   2.883  0.00556 ** 
## thick_r      0.92174    0.08979 82.77441  10.266  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## thick_r -0.656

Volume

volume_reg_ir <- lmer(vol_i ~ vol_r * species + (1 | individual_id), data = volume_all)
summary(volume_reg_ir)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_r * species + (1 | individual_id)
##    Data: volume_all
## 
## REML criterion at convergence: 680.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.1113 -0.2169 -0.0116  0.2568  4.1150 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 1.765    1.329   
##  Residual                  3.864    1.966   
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         1.30390    0.88276  67.59939   1.477    0.144    
## vol_r               1.43146    0.06500 137.46633  22.022  < 2e-16 ***
## speciesCC           0.09332    1.30665  68.25035   0.071    0.943    
## speciesCO          -1.05684    1.30904  58.95997  -0.807    0.423    
## speciesCYOS         0.72007    1.50989  99.54528   0.477    0.634    
## speciesDL           4.19610    5.94891 121.58941   0.705    0.482    
## speciesDP          -0.04726    1.13820  72.30213  -0.042    0.967    
## speciesEGME        -0.97200    1.28596  66.89426  -0.756    0.452    
## speciesGR          -3.05725    2.23710  70.34650  -1.367    0.176    
## speciesNAND        -1.15642    2.09422  90.80367  -0.552    0.582    
## speciesPOLA        -1.11325    2.86861 122.30784  -0.388    0.699    
## speciesR           -0.78440    1.38258  68.02105  -0.567    0.572    
## vol_r:speciesCC    -0.52246    0.09625 138.33744  -5.428 2.48e-07 ***
## vol_r:speciesCO    -0.69817    1.14178  81.96921  -0.611    0.543    
## vol_r:speciesCYOS   0.17432    0.51313 120.21491   0.340    0.735    
## vol_r:speciesDL    -1.93146    2.40841 100.95555  -0.802    0.424    
## vol_r:speciesDP     0.37187    0.28946 139.31565   1.285    0.201    
## vol_r:speciesEGME  -0.55641    0.79461 123.35933  -0.700    0.485    
## vol_r:speciesGR    -0.02982    0.20193 105.54459  -0.148    0.883    
## vol_r:speciesNAND   0.69625    1.54659 139.62288   0.450    0.653    
## vol_r:speciesPOLA  -0.02907    1.19204 101.05446  -0.024    0.981    
## vol_r:speciesR     -0.28390    0.51591 107.55214  -0.550    0.583    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
volume_reg_ir_nospecies <- lmer(vol_i ~ vol_r + (1 | individual_id), data = volume_all)
summary(volume_reg_ir_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_r + (1 | individual_id)
##    Data: volume_all
## 
## REML criterion at convergence: 770.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9613 -0.3012 -0.1227  0.2502  5.3131 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 2.812    1.677   
##  Residual                  4.812    2.194   
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   1.15637    0.32749  65.60806   3.531 0.000763 ***
## vol_r         1.21440    0.04462 155.39471  27.216  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## vol_r -0.459

b. Initial values predicted by herbarium values

This section includes code to fit linear mixed effect models to answer the question: how well do herbarium trait values predict initial trait values? Each trait is predicted by two models: initial values as a function of herbarium values with a species interaction and intiial values as a function of herbarium values without a species interaction - both models have individual_id included as a random effect. We view model summaries using summary().

Height

height_reg_ih <- lmer(height_i ~ height_h * species + (1 | individual_id), data = height_all)
summary(height_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: height_i ~ height_h * species + (1 | individual_id)
##    Data: height_all
## 
## REML criterion at convergence: 752
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3839 -0.3069 -0.0696  0.1275  7.4336 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.04431  0.2105  
##  Residual                  5.53471  2.3526  
## Number of obs: 167, groups:  individual_id, 54
## 
## Fixed effects:
##                       Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            0.53223    1.52878 101.18822   0.348   0.7285    
## height_h               1.08408    0.10738 117.09633  10.095   <2e-16 ***
## speciesCC             -2.62401    2.44318 134.06144  -1.074   0.2847    
## speciesCO             -0.61692    2.23952  74.89845  -0.275   0.7837    
## speciesCYOS            1.36457    3.52820 137.61217   0.387   0.6995    
## speciesDL              5.36919    3.03593 144.86101   1.769   0.0791 .  
## speciesDP              0.78055    1.92977  96.65881   0.404   0.6868    
## speciesEGME           -0.31909    1.83980  96.77177  -0.173   0.8627    
## speciesGR              1.50188    5.33634 112.08412   0.281   0.7789    
## speciesNAND            0.04037    2.76399 141.33926   0.015   0.9884    
## speciesPOLA           -0.15879    3.19891 143.87363  -0.050   0.9605    
## speciesR              -0.20268    2.04852  91.93949  -0.099   0.9214    
## height_h:speciesCC     0.25621    0.15139 133.99473   1.692   0.0929 .  
## height_h:speciesCO    -0.01880    0.21540  72.49721  -0.087   0.9307    
## height_h:speciesCYOS   0.05317    0.16504 136.64447   0.322   0.7478    
## height_h:speciesDL    -0.19227    0.15126 144.71505  -1.271   0.2057    
## height_h:speciesDP     0.04072    0.13133 110.99089   0.310   0.7571    
## height_h:speciesEGME  -0.01431    0.12181 126.89901  -0.117   0.9067    
## height_h:speciesGR    -0.11581    0.25850 104.85633  -0.448   0.6551    
## height_h:speciesNAND  -0.03487    0.24531 135.14081  -0.142   0.8872    
## height_h:speciesPOLA   0.02274    0.25690 120.08051   0.089   0.9296    
## height_h:speciesR      0.01836    0.16646  93.30354   0.110   0.9124    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#This model describes height_i as a function of height_h, with a species interaction term and random effect of individual_id

height_reg_ih_nospecies <- lmer(height_i ~ height_h + (1 | individual_id), data = height_all)
summary(height_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: height_i ~ height_h + (1 | individual_id)
##    Data: height_all
## 
## REML criterion at convergence: 783.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.9741 -0.4023 -0.1279  0.1446  7.3759 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.8207   0.9059  
##  Residual                  5.5105   2.3474  
## Number of obs: 167, groups:  individual_id, 54
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   0.25808    0.46707 115.71914   0.553    0.582    
## height_h      1.13738    0.02731 139.85308  41.641   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## height_h -0.880
#This model describes height_i as a function of height_h with a random effect of individual_id and without a species interaction term

Width

width_reg_ih <- lmer(width_i ~ width_h * species + (1 | individual_id), data = width_all)
summary(width_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_h * species + (1 | individual_id)
##    Data: width_all
## 
## REML criterion at convergence: 220.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4570 -0.3541 -0.0043  0.4800  3.7396 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.02578  0.1606  
##  Residual                  0.17125  0.4138  
## Number of obs: 162, groups:  individual_id, 54
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           0.272115   0.236478 126.740300   1.151 0.252021    
## width_h               1.151373   0.083930 138.419124  13.718  < 2e-16 ***
## speciesCC            -0.620575   0.342832 123.037547  -1.810 0.072714 .  
## speciesCO             0.446118   0.415375 139.326962   1.074 0.284674    
## speciesCYOS           0.471786   0.340900 109.731366   1.384 0.169185    
## speciesDL             0.246742   1.344039 112.207950   0.184 0.854673    
## speciesDP            -0.523920   0.296015 118.418792  -1.770 0.079316 .  
## speciesEGME           0.284999   0.408095 132.278112   0.698 0.486176    
## speciesGR            -0.122673   0.502592 136.144809  -0.244 0.807536    
## speciesNAND          -0.126218   0.356086 102.952743  -0.354 0.723719    
## speciesPOLA          -0.433035   1.673731 107.929898  -0.259 0.796341    
## speciesR             -0.123335   0.301637 110.557970  -0.409 0.683415    
## width_h:speciesCC     0.182215   0.088658 138.359422   2.055 0.041735 *  
## width_h:speciesCO    -0.288921   0.122342 134.718497  -2.362 0.019630 *  
## width_h:speciesCYOS  -0.090180   0.097947 139.835801  -0.921 0.358791    
## width_h:speciesDL    -0.373692   1.575169 100.246932  -0.237 0.812956    
## width_h:speciesDP     0.467741   0.135512 135.940672   3.452 0.000743 ***
## width_h:speciesEGME  -0.196477   0.145111 139.328362  -1.354 0.177934    
## width_h:speciesGR     0.008776   0.243724 134.944567   0.036 0.971328    
## width_h:speciesNAND  -0.053593   0.107060 134.929958  -0.501 0.617479    
## width_h:speciesPOLA   0.158052   0.600841 101.038027   0.263 0.793046    
## width_h:speciesR     -0.028081   0.091096 139.271757  -0.308 0.758347    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
width_reg_ih_nospecies <- lmer(width_i ~ width_h + (1 | individual_id), data = width_all)
summary(width_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: width_i ~ width_h + (1 | individual_id)
##    Data: width_all
## 
## REML criterion at convergence: 281.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7952 -0.4539  0.0449  0.3985  4.8626 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.1262   0.3552  
##  Residual                  0.2309   0.4805  
## Number of obs: 162, groups:  individual_id, 54
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   0.08737    0.08547  98.88700   1.022    0.309    
## width_h       1.20652    0.01892 154.81287  63.760   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## width_h -0.686

Perimeter

perimeter_reg_ih <- lmer(per_i ~ per_h * species + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_h * species + (1 | individual_id)
##    Data: perimeter_all
## 
## REML criterion at convergence: 1361.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0987 -0.2817 -0.0422  0.1080  4.7679 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 225.6    15.02   
##  Residual                  501.2    22.39   
## Number of obs: 153, groups:  individual_id, 52
## 
## Fixed effects:
##                    Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)         2.59684   13.77863  88.26129   0.188 0.850942    
## per_h               1.19298    0.10995 125.45680  10.850  < 2e-16 ***
## speciesCC          -6.82156   19.70908 103.74382  -0.346 0.729959    
## speciesCO         -19.35573   23.91988 114.13986  -0.809 0.420089    
## speciesCYOS         6.64731   21.07977 109.76272   0.315 0.753102    
## speciesDP          34.33115   17.91760  88.87197   1.916 0.058572 .  
## speciesEGME        -2.02207   19.61561 105.50238  -0.103 0.918091    
## speciesGR           1.05566   22.86537  80.61552   0.046 0.963290    
## speciesNAND       -19.63704   47.84044  78.89181  -0.410 0.682574    
## speciesPOLA       147.31861   41.95239 132.97267   3.512 0.000609 ***
## speciesR           -1.44377   18.60890  89.93418  -0.078 0.938330    
## per_h:speciesCC     0.11400    0.22636 126.28471   0.504 0.615398    
## per_h:speciesCO     0.18647    0.24758 132.38146   0.753 0.452687    
## per_h:speciesCYOS  -0.11782    0.13820 132.81018  -0.853 0.395437    
## per_h:speciesDP     0.10835    0.13834 122.61000   0.783 0.434995    
## per_h:speciesEGME  -0.09324    0.29545 110.83387  -0.316 0.752910    
## per_h:speciesGR     0.03776    0.19031 116.97364   0.198 0.843049    
## per_h:speciesNAND   0.20992    0.35881  95.46443   0.585 0.559902    
## per_h:speciesPOLA  -1.03952    0.28146 106.94185  -3.693 0.000351 ***
## per_h:speciesR      0.09234    0.17917 132.89084   0.515 0.607163    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
perimeter_reg_ih_nospecies <- lmer(per_i ~ per_h + (1 | individual_id), data = perimeter_all)
summary(perimeter_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: per_i ~ per_h + (1 | individual_id)
##    Data: perimeter_all
## 
## REML criterion at convergence: 1472.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.8908 -0.3120 -0.0885  0.1870  4.6894 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 563.2    23.73   
##  Residual                  573.3    23.94   
## Number of obs: 153, groups:  individual_id, 52
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  10.01440    5.75140 103.30764   1.741   0.0846 .  
## per_h         1.21717    0.04487 150.79674  27.126   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## per_h -0.740

Surface Area

surf_reg_ih <- lmer(surf_i ~ surf_h * species + (1 | individual_id), data = surf_all)
summary(surf_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_h * species + (1 | individual_id)
##    Data: surf_all
## 
## REML criterion at convergence: 1146.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0489 -0.2322  0.0118  0.2319  3.3970 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 99.13    9.956   
##  Residual                  80.88    8.993   
## Number of obs: 155, groups:  individual_id, 52
## 
## Fixed effects:
##                      Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         -0.463120   7.163540 101.638932  -0.065  0.94858    
## surf_h               1.652726   0.110907 133.491909  14.902  < 2e-16 ***
## speciesCC            1.695839   9.223655  86.168702   0.184  0.85456    
## speciesCO            2.335557   9.411117  97.867238   0.248  0.80452    
## speciesCYOS          2.814195  10.432608 112.263576   0.270  0.78785    
## speciesDP           15.162740   8.453638  97.601006   1.794  0.07597 .  
## speciesEGME          1.351731   9.298801  88.467416   0.145  0.88475    
## speciesGR           -1.999914  14.572854  81.746029  -0.137  0.89118    
## speciesNAND          1.678109  14.062657 115.197954   0.119  0.90522    
## speciesPOLA        -87.262355  25.937553 132.709192  -3.364  0.00100 ** 
## speciesR             0.702535  10.027528  79.063684   0.070  0.94432    
## surf_h:speciesCC    -0.042041   0.115345 132.677458  -0.364  0.71608    
## surf_h:speciesCO    -0.390872   0.504958 112.882983  -0.774  0.44051    
## surf_h:speciesCYOS  -0.042381   0.160204 119.030005  -0.265  0.79182    
## surf_h:speciesDP     0.007247   0.125121 133.330151   0.058  0.95390    
## surf_h:speciesEGME  -0.435098   0.142602 123.731498  -3.051  0.00279 ** 
## surf_h:speciesGR    -0.212427   0.206861 113.434524  -1.027  0.30665    
## surf_h:speciesNAND  -0.064651   0.338893 108.339360  -0.191  0.84906    
## surf_h:speciesPOLA   2.188812   0.450745  99.174844   4.856 4.47e-06 ***
## surf_h:speciesR     -0.245841   0.197231 105.183512  -1.246  0.21536    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
surf_reg_ih_nospecies <- lmer(surf_i ~ surf_h + (1 | individual_id), data = surf_all)
summary(surf_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: surf_i ~ surf_h + (1 | individual_id)
##    Data: surf_all
## 
## REML criterion at convergence: 1261.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8828 -0.3740 -0.0307  0.3144  3.5874 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 187.6    13.70   
##  Residual                  111.4    10.56   
## Number of obs: 155, groups:  individual_id, 52
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)   2.33827    2.47842  81.86394   0.943    0.348    
## surf_h        1.59703    0.02794 140.02936  57.160   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr)
## surf_h -0.532

Weight

tdmc_reg_ih <- lmer(weight_i ~ weight_h * species + (1 | individual_id), data = tdmc_all)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(tdmc_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_h * species + (1 | individual_id)
##    Data: tdmc_all
## 
## REML criterion at convergence: 2630.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9009 -0.1250 -0.0089  0.1327  7.8408 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept)  466323   682.9  
##  Residual                  2460787  1568.7  
## Number of obs: 159, groups:  individual_id, 55
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -112.6370   598.8637    71.1358  -0.188 0.851346    
## weight_h                 7.8059     0.2498   138.8975  31.246  < 2e-16 ***
## speciesCC              205.9332   927.2423    77.4660   0.222 0.824827    
## speciesCO              -23.2533   954.1125    74.7643  -0.024 0.980621    
## speciesCYOS            227.2907  1276.3765   130.1371   0.178 0.858941    
## speciesDP              155.9781   798.6326    70.8440   0.195 0.845712    
## speciesEGME            168.0942   940.4295    83.5829   0.179 0.858573    
## speciesGR            -1545.9601  1472.0010    68.2519  -1.050 0.297311    
## speciesNAND           2311.3346  1718.3377   134.9808   1.345 0.180849    
## speciesPOLA            129.8183  1807.9914   117.0949   0.072 0.942882    
## speciesR               223.2403   914.6747    80.4121   0.244 0.807802    
## weight_h:speciesCC      -3.1801     0.3491   136.7897  -9.110 9.01e-16 ***
## weight_h:speciesCO      -5.3849     1.5817   117.6456  -3.405 0.000907 ***
## weight_h:speciesCYOS    -1.3612     1.2933   138.7170  -1.053 0.294373    
## weight_h:speciesDP      -1.0085     0.7308   112.6685  -1.380 0.170294    
## weight_h:speciesEGME    -0.6318     3.8956   133.7570  -0.162 0.871404    
## weight_h:speciesGR      -3.6201     0.3985   121.2700  -9.085 2.39e-15 ***
## weight_h:speciesNAND    -8.2313     3.9984   115.6620  -2.059 0.041773 *  
## weight_h:speciesPOLA    -0.1031     4.1748   100.1666  -0.025 0.980344    
## weight_h:speciesR       -3.4064     1.1599    96.2291  -2.937 0.004151 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling
tdmc_reg_ih_nospecies <- lmer(weight_i ~ weight_h + (1 | individual_id), data = tdmc_all)
## Warning: Some predictor variables are on very different scales: consider
## rescaling
## Warning: Some predictor variables are on very different scales: consider
## rescaling
summary(tdmc_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: weight_i ~ weight_h + (1 | individual_id)
##    Data: tdmc_all
## 
## REML criterion at convergence: 2959.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8499 -0.1521 -0.0275  0.1360  7.5307 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 4023996  2006    
##  Residual                  5213112  2283    
## Number of obs: 159, groups:  individual_id, 55
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept) 157.2287   382.9288  47.7765   0.411    0.683    
## weight_h      5.2341     0.1939 139.7789  26.997   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## weight_h -0.504
## fit warnings:
## Some predictor variables are on very different scales: consider rescaling

Thickness

thick_reg_ih <- lmer(thick_i ~ thick_h * species + (1 | individual_id), data = thickness_all)
summary(thick_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_h * species + (1 | individual_id)
##    Data: thickness_all
## 
## REML criterion at convergence: -112.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.5330 -0.2257 -0.0360  0.1393  7.2536 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.04748  0.2179  
##  Residual                  0.01649  0.1284  
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##                      Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)           0.11225    0.11567  76.23562   0.970   0.3349    
## thick_h               0.47944    0.86955 107.00347   0.551   0.5825    
## speciesCC             0.31780    0.16174  63.87740   1.965   0.0538 .  
## speciesCO             0.31043    0.35703 110.60168   0.869   0.3865    
## speciesCYOS           0.16492    0.30022 139.77289   0.549   0.5837    
## speciesDL             0.05775    0.73919 121.33255   0.078   0.9379    
## speciesDP             0.50291    0.22373 130.49152   2.248   0.0263 *  
## speciesEGME           0.03790    0.18241  90.42927   0.208   0.8359    
## speciesGR             3.77862    0.82284 139.99471   4.592 9.68e-06 ***
## speciesNAND          -0.01503    0.92112 134.73544  -0.016   0.9870    
## speciesPOLA          -0.06745    0.38776 128.17496  -0.174   0.8622    
## speciesR             -0.10166    0.20280  60.71286  -0.501   0.6180    
## thick_h:speciesCC    -0.31147    0.90024 106.65416  -0.346   0.7300    
## thick_h:speciesCO    -0.24292    1.08623 133.32579  -0.224   0.8234    
## thick_h:speciesCYOS   1.39921    1.55498 109.97012   0.900   0.3702    
## thick_h:speciesDL    -0.47944    5.44469  98.06901  -0.088   0.9300    
## thick_h:speciesDP    -2.79756    1.73680 129.77461  -1.611   0.1097    
## thick_h:speciesEGME   0.13574    2.48680 105.53951   0.055   0.9566    
## thick_h:speciesGR    -2.41704    1.10667 123.25325  -2.184   0.0309 *  
## thick_h:speciesNAND   0.28492   11.28077 137.55240   0.025   0.9799    
## thick_h:speciesPOLA   1.40537    6.71613  97.98956   0.209   0.8347    
## thick_h:speciesR      1.21838    1.43063 115.64928   0.852   0.3962    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
thick_reg_ih_nospecies <- lmer(thick_i ~ thick_h + (1 | individual_id), data = thickness_all)
summary(thick_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: thick_i ~ thick_h + (1 | individual_id)
##    Data: thickness_all
## 
## REML criterion at convergence: -47.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1179 -0.2023 -0.0820  0.1024  6.7408 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.05453  0.2335  
##  Residual                  0.02079  0.1442  
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.20382    0.04228 47.93657   4.821 1.48e-05 ***
## thick_h      0.95958    0.09883 75.47472   9.710 6.32e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##         (Intr)
## thick_h -0.588

Volume

volume_reg_ih <- lmer(vol_i ~ vol_h * species + (1 | individual_id), data = volume_all)
summary(volume_reg_ih)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_h * species + (1 | individual_id)
##    Data: volume_all
## 
## REML criterion at convergence: 722.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8208 -0.2484 -0.0213  0.2188  4.4006 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 0.3959   0.6292  
##  Residual                  7.2424   2.6912  
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##                   Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)        -0.2386     1.0413  97.6340  -0.229 0.819210    
## vol_h               6.1714     0.3727 139.9985  16.557  < 2e-16 ***
## speciesCC           0.1428     1.5687  96.6246   0.091 0.927664    
## speciesCO           0.4539     1.4577  87.4885   0.311 0.756233    
## speciesCYOS         0.8583     2.2830 139.9888   0.376 0.707525    
## speciesDL           5.2386     4.8174 131.1702   1.087 0.278835    
## speciesDP           0.9271     1.3266  87.0716   0.699 0.486484    
## speciesEGME         0.7993     1.4005  79.3197   0.571 0.569798    
## speciesGR          -1.5082     2.7201 107.8951  -0.554 0.580414    
## speciesNAND         1.2685     2.5855 133.2270   0.491 0.624496    
## speciesPOLA         3.7132     3.7907 139.3476   0.980 0.328996    
## speciesR            0.3816     1.7283  81.9842   0.221 0.825814    
## vol_h:speciesCC    -1.5920     0.5894 128.1774  -2.701 0.007845 ** 
## vol_h:speciesCO    -5.3670     1.3987 124.1825  -3.837 0.000197 ***
## vol_h:speciesCYOS  -2.1893     1.7567 117.9363  -1.246 0.215151    
## vol_h:speciesDL    -6.6714     3.3170 106.3239  -2.011 0.046827 *  
## vol_h:speciesDP    -2.0224     0.8121 108.2204  -2.490 0.014280 *  
## vol_h:speciesEGME  -4.9200     2.0074 136.3490  -2.451 0.015517 *  
## vol_h:speciesGR    -3.6149     0.5987 137.4855  -6.038 1.38e-08 ***
## vol_h:speciesNAND  -3.9645     4.1135 127.1705  -0.964 0.336993    
## vol_h:speciesPOLA  -7.6968     7.0171 105.8129  -1.097 0.275191    
## vol_h:speciesR     -2.9773     1.8410  68.9227  -1.617 0.110406    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
volume_reg_ih_nospecies <- lmer(vol_i ~ vol_h + (1 | individual_id), data = volume_all)
summary(volume_reg_ih_nospecies)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: vol_i ~ vol_h + (1 | individual_id)
##    Data: volume_all
## 
## REML criterion at convergence: 888
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0468 -0.2718 -0.0058  0.2196  5.9598 
## 
## Random effects:
##  Groups        Name        Variance Std.Dev.
##  individual_id (Intercept) 6.823    2.612   
##  Residual                  9.825    3.135   
## Number of obs: 162, groups:  individual_id, 53
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)   0.2726     0.5367  57.1506   0.508    0.613    
## vol_h         3.7867     0.2330 142.9514  16.252   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##       (Intr)
## vol_h -0.573

#3. Comparing models using AIC

This section includes code comparing the AIC values for models with a species interaction term and without a species interaction term to answer the question: Do models with a species interaction improve the model’s fit? All trait models for rehydrated and herbarium traits are tested, and the results are commented below each model.

a. Comparisons for models with rehydrated traits as predictors

Height

#..........................Comparisons for Height Models........................

model.sel(height_reg_ir, height_reg_ir_nospecies)
## Model selection table 
##                            (Int) hgh_r spc hgh_r:spc df   logLik  AICc delta
## height_reg_ir_nospecies  0.04065 1.024                4 -169.859 348.0  0.00
## height_reg_ir           -0.27440 1.029   +         + 24 -169.996 396.6 48.66
##                         weight
## height_reg_ir_nospecies      1
## height_reg_ir                0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#nospecies is better

Width

#..........................Comparisons for Width Models.........................

model.sel(width_reg_ir, width_reg_ir_nospecies)
## Model selection table 
##                          (Int) spc  wdt_r spc:wdt_r df logLik  AICc  delta
## width_reg_ir_nospecies 0.08645     0.9939            4 19.245 -30.2   0.00
## width_reg_ir           0.00229   + 1.0300         + 24 -8.349  73.5 103.69
##                        weight
## width_reg_ir_nospecies      1
## width_reg_ir                0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#nospecies is better

Perimeter

#..........................Comparisons for Perimeter Models......................

model.sel(perimeter_reg_ir, perimeter_reg_ir_nospecies)
## Model selection table 
##                             (Int) per_r spc per_r:spc df   logLik   AICc delta
## perimeter_reg_ir           11.940 0.899   +         + 22 -649.166 1350.1  0.00
## perimeter_reg_ir_nospecies  5.569 1.007                4 -717.235 1442.7 92.62
##                            weight
## perimeter_reg_ir                1
## perimeter_reg_ir_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#with species is better

Surface Area

#..........................Comparisons for Surface Area Models..................

model.sel(surf_reg_ir, surf_reg_ir_nospecies)
## Model selection table 
##                        (Int) spc srf_r spc:srf_r df   logLik   AICc delta
## surf_reg_ir           -1.581   + 1.129         + 22 -594.245 1240.2  0.00
## surf_reg_ir_nospecies  8.157     1.038            4 -636.504 1281.3 41.12
##                       weight
## surf_reg_ir                1
## surf_reg_ir_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#species is better

Weight

#..........................Comparisons for Weight Models........................

model.sel(tdmc_reg_ir, tdmc_reg_ir_nospecies)
## Model selection table 
##                        (Int) spc wgh_r spc:wgh_r df    logLik   AICc  delta
## tdmc_reg_ir           -859.7   + 1.809         + 22 -1278.379 2608.2   0.00
## tdmc_reg_ir_nospecies  832.5     1.398            4 -1453.517 2915.3 307.09
##                       weight
## tdmc_reg_ir                1
## tdmc_reg_ir_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#species is better

Thickness

#..........................Comparisons for Thickness Models.....................

model.sel(thick_reg_ir, thick_reg_ir_nospecies)
## Model selection table 
##                         (Int) spc  thc_r spc:thc_r df logLik  AICc delta weight
## thick_reg_ir_nospecies 0.1342     0.9217            4 31.946 -55.6  0.00      1
## thick_reg_ir           0.1205   + 0.1399         + 24 48.267 -39.8 15.86      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#nospecies is better, but rehydrated thickness it's not a good predictor overall

Volume

#..........................Comparisons for Volume Models........................

model.sel(volume_reg_ir, volume_reg_ir_nospecies)
## Model selection table 
##                         (Int) spc vol_r spc:vol_r df   logLik  AICc delta
## volume_reg_ir           1.304   + 1.431         + 24 -340.348 737.5  0.00
## volume_reg_ir_nospecies 1.156     1.214            4 -385.284 778.8 41.37
##                         weight
## volume_reg_ir                1
## volume_reg_ir_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#with species is better

##b. Comparisons for models with Herbarium traits as a predictor variable

Height

#..........................Comparisons for Height Models........................

model.sel(height_reg_ih, height_reg_ih_nospecies)
## Model selection table 
##                          (Int) hgh_h spc hgh_h:spc df   logLik  AICc delta
## height_reg_ih_nospecies 0.2581 1.137                4 -391.826 791.9  0.00
## height_reg_ih           0.5322 1.084   +         + 24 -376.005 808.5 16.56
##                         weight
## height_reg_ih_nospecies      1
## height_reg_ih                0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#nospecies is better

Width

#..........................Comparisons for Width Models........................

model.sel(width_reg_ih, width_reg_ih_nospecies)
## Model selection table 
##                          (Int) spc wdt_h spc:wdt_h df   logLik  AICc delta
## width_reg_ih           0.27210   + 1.151         + 24 -110.139 277.0   0.0
## width_reg_ih_nospecies 0.08737     1.207            4 -140.689 289.6  12.6
##                        weight
## width_reg_ih            0.998
## width_reg_ih_nospecies  0.002
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#species is better

Perimeter

#.......................Comparisons for Perimeter Models........................

model.sel(perimeter_reg_ih, perimeter_reg_ih_nospecies)
## Model selection table 
##                             (Int) per_h spc per_h:spc df   logLik   AICc delta
## perimeter_reg_ih            2.597 1.193   +         + 22 -680.895 1413.6  0.00
## perimeter_reg_ih_nospecies 10.010 1.217                4 -736.285 1480.8 67.27
##                            weight
## perimeter_reg_ih                1
## perimeter_reg_ih_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#with species is better

Surface Area

#....................Comparisons for Surface Area Models........................

model.sel(surf_reg_ih, surf_reg_ih_nospecies)
## Model selection table 
##                         (Int) spc srf_h spc:srf_h df   logLik   AICc delta
## surf_reg_ih           -0.4631   + 1.653         + 22 -573.459 1198.6  0.00
## surf_reg_ih_nospecies  2.3380     1.597            4 -630.864 1270.0 71.41
##                       weight
## surf_reg_ih                1
## surf_reg_ih_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#species is better

Weight

#..........................Comparisons for Weight Models........................

model.sel(tdmc_reg_ih, tdmc_reg_ih_nospecies)
## Model selection table 
##                        (Int) spc wgh_h spc:wgh_h df    logLik   AICc  delta
## tdmc_reg_ih           -112.6   + 7.806         + 22 -1315.397 2682.2   0.00
## tdmc_reg_ih_nospecies  157.2     5.234            4 -1479.581 2967.4 285.19
##                       weight
## tdmc_reg_ih                1
## tdmc_reg_ih_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#species is better

Thickness

#.......................Comparisons for Thickness Models........................

model.sel(thick_reg_ih, thick_reg_ih_nospecies)
## Model selection table 
##                         (Int) spc  thc_h spc:thc_h df logLik  AICc delta weight
## thick_reg_ih           0.1123   + 0.4794         + 24 56.121 -55.5  0.00      1
## thick_reg_ih_nospecies 0.2038     0.9596            4 23.776 -39.3 16.19      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#species is better, but doesn't really matter cause it's not a good predictor

Volume

#..........................Comparisons for Volume Models........................

model.sel(volume_reg_ih, volume_reg_ih_nospecies)
## Model selection table 
##                           (Int) spc vol_h spc:vol_h df   logLik  AICc  delta
## volume_reg_ih           -0.2386   + 6.171         + 24 -361.033 778.8   0.00
## volume_reg_ih_nospecies  0.2726     3.787            4 -443.999 896.3 117.43
##                         weight
## volume_reg_ih                1
## volume_reg_ih_nospecies      0
## Models ranked by AICc(x) 
## Random terms (all models): 
##   1 | individual_id
#with species is better

4. Calculating R squared values for the best fitting models

This section includes code that calculates the marginal and conditional R squared values of the best fitting models from the section above. We calculate the R squared value using MuMIn::r.squaredGLMM(model_name) from the MuMIn package.

##a. Calculating R squared values for models with rehydrated traits as a predictor

Height

#..................R^2 for Height Model without species interaction term........

MuMIn::r.squaredGLMM(height_reg_ir_nospecies)
##            R2m       R2c
## [1,] 0.9933504 0.9959516

Width

#..................R^2 for Width Model without species interaction term........

MuMIn::r.squaredGLMM(width_reg_ir_nospecies)
##            R2m       R2c
## [1,] 0.9960237 0.9960237

Perimeter

#...............R^2 for Perimeter Model with species interaction term........

MuMIn::r.squaredGLMM(perimeter_reg_ir)
##            R2m       R2c
## [1,] 0.9417458 0.9552199

Surface Area

#............R^2 for Surface Area Model with species interaction term........

MuMIn::r.squaredGLMM(surf_reg_ir)
##            R2m       R2c
## [1,] 0.9630645 0.9782143

Weight

#..................R^2 for Weight Model with species interaction term........

MuMIn::r.squaredGLMM(tdmc_reg_ir)
##            R2m       R2c
## [1,] 0.9728602 0.9728602

Thickness

#...............R^2 for Thickness Model without species interaction term........

MuMIn::r.squaredGLMM(thick_reg_ir_nospecies)
##            R2m       R2c
## [1,] 0.4802988 0.8816962

Volume

#..................R^2 for Volume Model with species interaction term........

MuMIn::r.squaredGLMM(volume_reg_ir)
##            R2m       R2c
## [1,] 0.8780037 0.9162553

##b. Calculating R squared values for models with Herbarium traits as predictor variables

Height

#..................R^2 for Height Model without species interaction term........

MuMIn::r.squaredGLMM(height_reg_ih_nospecies)
##            R2m       R2c
## [1,] 0.9182208 0.9288213

Width

#..................R^2 for Width Model with species interaction term............

MuMIn::r.squaredGLMM(width_reg_ih)
##            R2m       R2c
## [1,] 0.9814827 0.9839057

Perimeter

#...............R^2 for Perimeter Model with species interaction term........

MuMIn::r.squaredGLMM(perimeter_reg_ih)
##            R2m       R2c
## [1,] 0.8997358 0.9308613

Surface Area

#............R^2 for Surface Area Model with species interaction term........

MuMIn::r.squaredGLMM(surf_reg_ih)
##            R2m       R2c
## [1,] 0.9662619 0.9848414

Weight

#..................R^2 for Weight Model with species interaction term........

MuMIn::r.squaredGLMM(tdmc_reg_ih)
##            R2m       R2c
## [1,] 0.9435099 0.9525094

Thickness

#..................R^2 for Thickness Model with species interaction term........

MuMIn::r.squaredGLMM(thick_reg_ih)
##            R2m       R2c
## [1,] 0.6413815 0.9075627

Volume

#..................R^2 for Volume Model with species interaction term...........

MuMIn::r.squaredGLMM(volume_reg_ih)
##            R2m       R2c
## [1,] 0.8312203 0.8399674

#5. Saving models as predictions for later plotting

This section includes code to save the model predictions of the best fitting models as data frames using ggpredict(). This includes predictions for all traits and both rehydrate and herbarium treatments.

a. Saving predictions for models with rehydrated traits as predictor variables

Height

     #........Predictions for the Height Model Without Species........

height_ir_predictions <- ggpredict(height_reg_ir_nospecies, terms = "height_r[0:45]") |> 
  #saving predictions for height_r values from 0-45 
  mutate(trait = "height") #adding a column to label the type of trait values

height_ir_predictions
## # Predicted values of height_i
## 
## height_r | Predicted |      95% CI |  trait
## -------------------------------------------
##        0 |      0.04 | -0.23, 0.31 | height
##        1 |      1.06 |  0.81, 1.32 | height
##        2 |      2.09 |  1.84, 2.34 | height
##        3 |      3.11 |  2.87, 3.35 | height
##        4 |      4.14 |  3.91, 4.36 | height
##        5 |      5.16 |  4.94, 5.38 | height
##        6 |      6.18 |  5.98, 6.39 | height
## 
## Adjusted for:
## * individual_id = 0 (population-level)

Width

     #........Predictions for the Width Model Without Species........

width_ir_predictions <- ggpredict(width_reg_ir_nospecies, terms = "width_r[0:30]") |> #max width = 24
  mutate(trait = "width") 

# look at the data frame
width_ir_predictions
## # Predicted values of width_i
## 
## width_r | Predicted |     95% CI | trait
## ----------------------------------------
##       0 |      0.09 | 0.04, 0.14 | width
##       1 |      1.08 | 1.04, 1.12 | width
##       2 |      2.07 | 2.04, 2.11 | width
##       3 |      3.07 | 3.04, 3.10 | width
##       4 |      4.06 | 4.03, 4.09 | width
##       5 |      5.06 | 5.02, 5.09 | width
##       6 |      6.05 | 6.01, 6.09 | width
## 
## Adjusted for:
## * individual_id = 0 (population-level)

Perimeter

     #........Predictions for the Height Model With Species........
perimeter_ir_predictions <- ggpredict(perimeter_reg_ir, terms = "per_r[0:375]") |> 
  mutate(trait = "perimeter")

     #...........Incorporating Species as a Prediction Term...........
perimeter_ir_predictions_species <- ggpredict(perimeter_reg_ir, terms = c("per_r", "species")) |> 
  as_tibble() |> #converting to a tibble
  rename(species = group)
# look at the data frame
perimeter_ir_predictions
## # Predicted values of per_i
## 
## per_r | Predicted |       95% CI |     trait
## --------------------------------------------
##     0 |     11.94 | -7.30, 31.19 | perimeter
##     1 |     12.84 | -6.30, 31.99 | perimeter
##     2 |     13.74 | -5.30, 32.79 | perimeter
##     3 |     14.64 | -4.30, 33.59 | perimeter
##     4 |     15.54 | -3.31, 34.39 | perimeter
##     5 |     16.44 | -2.31, 35.19 | perimeter
##     6 |     17.34 | -1.31, 35.99 | perimeter
## 
## Adjusted for:
## *       species = BF
## * individual_id = 0 (population-level)

Surface Area

     #........Predictions for the Height Model With Species........

surf_ir_predictions <- ggpredict(surf_reg_ir, terms = "surf_r[0:430]") |> 
  mutate(trait = "surface area")

     #...........Incorporating Species as a Prediction Term...........
surf_ir_predictions_species <- ggpredict(surf_reg_ir, terms = c("surf_r", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
surf_ir_predictions
## # Predicted values of surf_i
## 
## surf_r | Predicted |        95% CI |        trait
## -------------------------------------------------
##      0 |     -1.58 | -17.03, 13.87 | surface area
##      1 |     -0.45 | -15.76, 14.86 | surface area
##      2 |      0.68 | -14.49, 15.85 | surface area
##      3 |      1.81 | -13.22, 16.84 | surface area
##      4 |      2.94 | -11.96, 17.83 | surface area
##      5 |      4.06 | -10.69, 18.82 | surface area
##      6 |      5.19 |  -9.42, 19.81 | surface area
## 
## Adjusted for:
## *       species = BF
## * individual_id = 0 (population-level)

Weight

     #........Predictions for the Weight Model With Species........

tdmc_ir_predictions <- ggpredict(tdmc_reg_ir, terms = "weight_r[0:29050]") |> 
  as_tibble() |> 
  mutate(trait = "thallus dry matter content")

     #...........Incorporating Species as a Prediction Term...........
tdmc_ir_predictions_species <- ggpredict(tdmc_reg_ir, terms = c("weight_r", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
tdmc_ir_predictions
## # A tibble: 29,051 × 7
##        x predicted std.error conf.low conf.high group trait                     
##    <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct> <chr>                     
##  1     0     -860.      390.   -1630.     -89.4 1     thallus dry matter content
##  2     1     -858.      390.   -1628.     -87.6 1     thallus dry matter content
##  3     2     -856.      390.   -1626.     -85.8 1     thallus dry matter content
##  4     3     -854.      389.   -1624.     -84.1 1     thallus dry matter content
##  5     4     -852.      389.   -1623.     -82.3 1     thallus dry matter content
##  6     5     -851.      389.   -1621.     -80.6 1     thallus dry matter content
##  7     6     -849.      389.   -1619.     -78.8 1     thallus dry matter content
##  8     7     -847.      389.   -1617.     -77.1 1     thallus dry matter content
##  9     8     -845.      389.   -1615.     -75.3 1     thallus dry matter content
## 10     9     -843.      389.   -1613.     -73.6 1     thallus dry matter content
## # ℹ 29,041 more rows

Thickness

     #........Predictions for the Thickness Model Without Species........

thick_ir_predictions <- ggpredict(thick_reg_ir_nospecies, terms = "thick_r[0:2]")  |> 
  mutate(trait = "thickness")

# look at the data frame
thick_ir_predictions
## # Predicted values of thick_i
## 
## : 1
## 
## thick_r | Predicted |     95% CI |     trait
## --------------------------------------------
##       0 |      0.13 | 0.04, 0.23 | thickness
##       1 |      1.06 | 0.92, 1.19 | thickness
##       2 |      1.98 | 1.68, 2.28 | thickness
## 
## Adjusted for:
## * individual_id = 0 (population-level)

Volume

     #........Predictions for the Volume Model With Species........

volume_ir_predictions <- ggpredict(volume_reg_ir, terms = "vol_r[0:30]") |> 
  mutate(trait = "volume")

     #...........Incorporating Species as a Prediction Term...........
volume_ir_predictions_species <- ggpredict(volume_reg_ir, terms = c("vol_r", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
volume_ir_predictions
## # Predicted values of vol_i
## 
## vol_r | Predicted |       95% CI |  trait
## -----------------------------------------
##     0 |      1.30 | -0.44,  3.05 | volume
##     1 |      2.74 |  1.06,  4.41 | volume
##     2 |      4.17 |  2.56,  5.77 | volume
##     3 |      5.60 |  4.05,  7.15 | volume
##     4 |      7.03 |  5.53,  8.53 | volume
##     5 |      8.46 |  7.00,  9.92 | volume
##     6 |      9.89 |  8.46, 11.33 | volume
## 
## Adjusted for:
## *       species = BF
## * individual_id = 0 (population-level)

##b. Saving predictions for models with herbarium traits as predictor variables

Height

     #........Predictions for the Height Model Without Species........
height_ih_predictions <- ggpredict(height_reg_ih_nospecies, terms = "height_h[0:45]") 

height_ih_predictions
## # Predicted values of height_i
## 
## height_h | Predicted |       95% CI
## -----------------------------------
##        0 |      0.26 | -0.66,  1.18
##        6 |      7.08 |  6.43,  7.74
##       11 |     12.77 | 12.28, 13.26
##       17 |     19.59 | 19.14, 20.04
##       23 |     26.42 | 25.80, 27.03
##       28 |     32.10 | 31.28, 32.93
##       34 |     38.93 | 37.82, 40.04
##       45 |     51.44 | 49.77, 53.11
## 
## Adjusted for:
## * individual_id = 0 (population-level)

Width

 #........Predictions for the Width Model Without Species........
width_ih_predictions <- ggpredict(width_reg_ih_nospecies, terms = "width_h[0:30]") #max width = 24

# look at the data frame
width_ih_predictions
## # Predicted values of width_i
## 
## width_h | Predicted |       95% CI
## ----------------------------------
##       0 |      0.09 | -0.08,  0.26
##       4 |      4.91 |  4.79,  5.04
##       7 |      8.53 |  8.34,  8.72
##      11 |     13.36 | 13.04, 13.68
##      15 |     18.19 | 17.72, 18.65
##      19 |     23.01 | 22.40, 23.62
##      23 |     27.84 | 27.08, 28.59
##      30 |     36.28 | 35.27, 37.30
## 
## Adjusted for:
## * individual_id = 0 (population-level)
     #...........Incorporating Species as a Prediction Term...........
width_ih_predictions_species <- ggpredict(width_reg_ih, terms = c("width_h", "species")) |> 
  as_tibble() |> 
  rename(species = group)

Perimeter

#........Predictions for the Perimeter Model With Species........
perimeter_ih_predictions <- ggpredict(perimeter_reg_ih, terms = "per_h[0:375]") 

     #...........Incorporating Species as a Prediction Term...........
perimeter_ih_predictions_species <- ggpredict(perimeter_reg_ih, terms = c("per_h", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
perimeter_ih_predictions_species
## # A tibble: 180 × 6
##        x predicted std.error conf.low conf.high species
##    <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>  
##  1     0     2.60       13.8    -24.7      29.9 BF     
##  2     0    -4.22       14.1    -32.1      23.7 CC     
##  3     0   -16.8        19.6    -55.4      21.9 CO     
##  4     0     9.24       16.0    -22.3      40.8 CYOS   
##  5     0    36.9        11.5     14.3      59.6 DP     
##  6     0     0.575      14.0    -27.0      28.2 EGME   
##  7     0     3.65       18.2    -32.4      39.8 GR     
##  8     0   -17.0        45.8   -108.       73.6 NAND   
##  9     0   150.         39.6     71.5     228.  POLA   
## 10     0     1.15       12.5    -23.6      25.9 R      
## # ℹ 170 more rows

Surface Area

 #........Predictions for the Surface Area Model With Species........
surf_ih_predictions <- ggpredict(surf_reg_ih, terms = "surf_h[0:430]") 

     #...........Incorporating Species as a Prediction Term...........
surf_ih_predictions_species <- ggpredict(surf_reg_ih, terms = c("surf_h", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
surf_ih_predictions_species
## # A tibble: 150 × 6
##        x predicted std.error conf.low conf.high species
##    <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>  
##  1     0    -0.463      7.16   -14.6       13.7 BF     
##  2     0     1.23       5.81   -10.3       12.7 CC     
##  3     0     1.87       6.10   -10.2       13.9 CO     
##  4     0     2.35       7.58   -12.7       17.4 CYOS   
##  5     0    14.7        4.49     5.82      23.6 DP     
##  6     0     0.889      5.93   -10.8       12.6 EGME   
##  7     0    -2.46      12.7    -27.6       22.6 GR     
##  8     0     1.21      12.1    -22.7       25.2 NAND   
##  9     0   -87.7       24.9   -137.       -38.4 POLA   
## 10     0     0.239      7.02   -13.6       14.1 R      
## # ℹ 140 more rows

Weight

     #........Predictions for the Weight Model With Species........

tdmc_ih_predictions <- ggpredict(tdmc_reg_ih, terms = "weight_h[0:29050]") 

     #...........Incorporating Species as a Prediction Term...........
tdmc_ih_predictions_species <- ggpredict(tdmc_reg_ih, terms = c("weight_h", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
tdmc_ih_predictions
## # Predicted values of weight_i
## 
## weight_h | Predicted |               95% CI
## -------------------------------------------
##        0 |   -112.64 |  -1296.85,   1071.57
##     3631 |  28230.75 |  26842.96,  29618.54
##     7263 |  56581.94 |  53600.77,  59563.12
##    10894 |  84925.33 |  80204.85,  89645.81
##    14525 | 113268.72 | 106779.23, 119758.20
##    18156 | 141612.10 | 133342.94, 149881.26
##    21787 | 169955.49 | 159901.65, 180009.33
##    29050 | 226650.07 | 213019.24, 240280.90
## 
## Adjusted for:
## *       species = BF
## * individual_id = 0 (population-level)
tdmc_ih_predictions_species
## # A tibble: 170 × 6
##        x predicted std.error conf.low conf.high species
##    <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>  
##  1     0    -113.       599.   -1297.     1072. BF     
##  2     0      93.3      708.   -1307.     1493. CC     
##  3     0    -136.       743.   -1605.     1333. CO     
##  4     0     115.      1127.   -2114.     2344. CYOS   
##  5     0      43.3      528.   -1001.     1088. DP     
##  6     0      55.5      725.   -1378.     1489. EGME   
##  7     0   -1659.      1345.   -4318.     1000. GR     
##  8     0    2199.      1611.    -986.     5384. NAND   
##  9     0      17.2     1706.   -3356.     3391. POLA   
## 10     0     111.       691.   -1257.     1478. R      
## # ℹ 160 more rows

Thickness

     #........Predictions for the Thickness Model Without Species........

thick_ih_predictions <- ggpredict(thick_reg_ih_nospecies, terms = "thick_h[0:2]") 

# look at the data frame
thick_ih_predictions
## # Predicted values of thick_i
## 
## thick_h | Predicted |     95% CI
## --------------------------------
##       0 |      0.20 | 0.12, 0.29
##       1 |      1.16 | 1.00, 1.32
##       2 |      2.12 | 1.78, 2.47
## 
## Adjusted for:
## * individual_id = 0 (population-level)
     #...........Incorporating Species as a Prediction Term...........
thick_ih_predictions_species <- ggpredict(thick_reg_ih, terms = c("thick_h", "species"))  |> 
  as_tibble() |> 
  rename(species = group)

Volume

     #........Predictions for the Volume Model With Species........
volume_ih_predictions <- ggpredict(volume_reg_ih, terms = "vol_h[0:30]") 

     #...........Incorporating Species as a Prediction Term...........
volume_ih_predictions_species <- ggpredict(volume_reg_ih, terms = c("vol_h", "species")) |> 
  as_tibble() |> 
  rename(species = group)
# look at the data frame
volume_ih_predictions
## # Predicted values of vol_i
## 
## vol_h | Predicted |         95% CI
## ----------------------------------
##     0 |     -0.24 |  -2.30,   1.82
##     4 |     24.45 |  22.50,  26.39
##     7 |     42.96 |  39.11,  46.81
##    11 |     67.65 |  60.95,  74.34
##    15 |     92.33 |  82.73, 101.93
##    19 |    117.02 | 104.49, 129.54
##    23 |    141.70 | 126.24, 157.16
##    30 |    184.90 | 164.30, 205.51
## 
## Adjusted for:
## *       species = BF
## * individual_id = 0 (population-level)
volume_ih_predictions_species
## # A tibble: 209 × 6
##        x predicted std.error conf.low conf.high species
##    <dbl>     <dbl>     <dbl>    <dbl>     <dbl> <fct>  
##  1   0.1     0.378     1.01    -1.63       2.38 BF     
##  2   0.1     0.362     1.14    -1.89       2.61 CC     
##  3   0.1     0.296     0.925   -1.53       2.12 CO     
##  4   0.1     1.02      1.87    -2.69       4.72 CYOS   
##  5   0.1     4.95      4.40    -3.74      13.6  DL     
##  6   0.1     1.10      0.765   -0.410      2.62 DP     
##  7   0.1     0.686     0.799   -0.894      2.27 EGME   
##  8   0.1    -1.49      2.47    -6.38       3.39 GR     
##  9   0.1     1.25      2.02    -2.75       5.25 NAND   
## 10   0.1     3.32      3.02    -2.65       9.29 POLA   
## # ℹ 199 more rows

6. Visualizing Models

This section includes code for scatterplots of model predictions using ggplot. Every model has one scatterplot with initial trait values on the Y axis and treatment (rehydrated/herbarium) trait values on the X axis. Models that incorporate a species interaction have an additional scatterplot faceted by species.

##a. Visualizing models with Rehydrated traits on the X axis

Height

     #.....Scatterplot of height_r(x) as Predictor of height_i(y).....

height_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = height_all,
             aes(x = height_r,
                 y = height_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = height_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = height_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate height (mm)",
       y = "Initial height (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 50, by = 10), 
                     limits = c(0, 50)) +
  scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10), 
                     limits = c(0, 50))+
  #add equation + r squared
   annotate("text", x = 10, y = 45, label= "y = 1.02x + 0.0407
R squared = 0.996")
  

height_ir_plot
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

### Width

     #.....Scatterplot of width_r(x) as Predictor of width_i(y).....

width_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = width_all,
             aes(x = width_r,
                 y = width_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = width_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = width_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate Width (mm)",
       y = "Initial width (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 30, by = 5), 
                     limits = c(0, 30)) +
  scale_y_continuous(breaks = seq(from = 0, to = 30, by = 5), 
                     limits = c(0, 30))+
  #add equation + r squared
   annotate("text", x = 10, y = 25, label= "y = 0.994x + 0.0865
R squared = 0.996")
  
width_ir_plot

### Thickness

     #.....Scatterplot of thickness_r(x) as Predictor of thickness_i(y).....

thick_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = thickness_all,
             aes(x = thick_r,
                 y = thick_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = thick_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = thick_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate thickness (mm)",
       y = "Initial thickness (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 2, by = 0.25), 
                     limits = c(0, 2)) +
  scale_y_continuous(breaks = seq(from = 0, to = 2, by = 0.25), 
                     limits = c(0, 2))+
  #add equation + r squared
   annotate("text", x = 0.75, y = 1.75, label= "y = 0.922x + 0.134
R squared = 0.882")

thick_ir_plot
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_line()`).

Perimeter

     #.....Scatterplot of perimeter_r(x) as Predictor of perimeter_i(y).....
perimeter_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = perimeter_all,
             aes(x = per_r,
                 y = per_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = perimeter_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = perimeter_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate perimeter (mm)",
       y = "Initial perimeter (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 375, by = 50), 
                     limits = c(0, 375)) +
  scale_y_continuous(breaks = seq(from = 0, to = 375, by = 50), 
                     limits = c(0, 375))+
  #add equation + r squared
   annotate("text", x = 75, y = 300, label= "y = 1.02x + 0.0407
R squared = 0.996")
perimeter_ir_plot
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

  #.Scatterplot of height_r(x) as Predictor of height_i(y) with a species facet.

perimeter_ir_plot_species <- ggplot() +
  
  geom_point(data = perimeter_all,
             aes(x = per_r,
                 y = per_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = perimeter_ir_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = perimeter_ir_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Rehydrate perimeter (mm)",
       y = "Initial perimeter (mm)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  perimeter_ir_plot_species
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

### Surface Area

     #.....Scatterplot of surf_r(x) as Predictor of surf_i(y).....

surf_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = surf_all,
             aes(x = surf_r,
                 y = surf_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = surf_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = surf_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate Surface Area (mm^2)",
       y = "Initial Surface Area (mm^2)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 450, by = 75), 
                     limits = c(0, 450)) +
  scale_y_continuous(breaks = seq(from = 0, to = 450, by = 75), 
                     limits = c(0, 450))
surf_ir_plot
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

 #.Scatterplot of surf_r(x) as Predictor of surf_i(y) with a species facet.
surf_ir_plot_species <- ggplot() +
  
  geom_point(data = surf_all,
             aes(x = surf_r,
                 y = surf_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = surf_ir_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = surf_ir_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Rehydrate surface area (mm^2)",
       y = "Initial surface area (mm^2)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  surf_ir_plot_species

### Weight

 #.....Scatterplot of weight_r(x) as Predictor of weight_i(y).....

tdmc_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = tdmc_all,
             aes(x = weight_r,
                 y = weight_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = tdmc_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = tdmc_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate weight (g)",
       y = "Initial weight (g)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 4850, by = 810), 
                     limits = c(0, 4850)) +
  scale_y_continuous(breaks = seq(from = 0, to = 4850, by = 810), 
                     limits = c(0, 4850))
tdmc_ir_plot
## Warning: Removed 54 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 26370 rows containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of weight_r(x) as Predictor of weight_i(y) with a species facet.


tdmc_ir_plot_species <- ggplot() +
  
  geom_point(data = tdmc_all,
             aes(x = weight_r,
                 y = weight_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = tdmc_ir_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = tdmc_ir_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Rehydrate weight (g)",
       y = "Initial weight (g)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


 tdmc_ir_plot_species
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

### Volume

#.....Scatterplot of vol_r(x) as Predictor of vol_i(y).....

volume_ir_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = volume_all,
             aes(x = vol_r,
                 y = vol_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = volume_ir_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = volume_ir_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate volume (mL)",
       y = "Initial volume (mL)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 15, by = 2.5), 
                     limits = c(0, 15)) +
  scale_y_continuous(breaks = seq(from = 0, to = 15, by = 2.5), 
                     limits = c(0, 15))
volume_ir_plot
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).

  #.Scatterplot of volume_r(x) as Predictor of volume_i(y) with a species facet.


volume_ir_plot_species <- ggplot() +
  
  geom_point(data = volume_all,
             aes(x = vol_r,
                 y = vol_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = volume_ir_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = volume_ir_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Rehydrate volume (mL)",
       y = "Initial volume (mL)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  volume_ir_plot_species

##b. Visualizing models with herbarium traits on the X axis

Height

     #....Scatterplot of height_h(x) as a Predictor of height_i(y)....
height_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = height_all,
             aes(x = height_h,
                 y = height_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = height_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = height_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Herbarium height (mm)",
       y = "Initial height (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 50, by = 10), 
                     limits = c(0, 50)) +
  scale_y_continuous(breaks = seq(from = 0, to = 50, by = 10), 
                     limits = c(0, 50))+
  #add equation + r squared
   annotate("text", x = 10, y = 45, label= "y = 1.14x + 0.258
R squared = 0.929")

height_ih_plot
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

Width

 #....Scatterplot of width_h(x) as a Predictor of width_i(y)....
width_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = width_all,
             aes(x = width_h,
                 y = width_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = width_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = width_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Herbarium width (mm)",
       y = "Initial width (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 20, by = 5), 
                     limits = c(0, 20)) +
  scale_y_continuous(breaks = seq(from = 0, to = 20, by = 5), 
                     limits = c(0, 20))
width_ih_plot
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of width_h(x) as Predictor of width_i(y) with a species facet.
width_ih_plot_species <- ggplot() +
  
  geom_point(data = width_all,
             aes(x = width_h,
                 y = width_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = width_ih_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = width_ih_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Herbarium width (mm)",
       y = "Initial width (mm)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  width_ih_plot_species

### Thickness

     #....Scatterplot of thick_h(x) as a Predictor of thick_i(y)....
thick_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = thickness_all,
             aes(x = thick_h,
                 y = thick_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = thick_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = thick_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Rehydrate thickness (mm)",
       y = "Initial thickness (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 2, by = 0.25), 
                     limits = c(0, 2)) +
  scale_y_continuous(breaks = seq(from = 0, to = 2, by = 0.25), 
                     limits = c(0, 2))

thick_ih_plot
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of thick_h(x) as Predictor of thick_i(y) with a species facet.
thick_ih_plot_species <- ggplot() +
  
  geom_point(data = thickness_all,
             aes(x = thick_h,
                 y = thick_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = thick_ih_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = thick_ih_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Herbarium thickness (mm)",
       y = "Initial thickness (mm)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  thick_ih_plot_species

### Perimeter

     #....Scatterplot of per_h(x) as a Predictor of per_i(y)....
perimeter_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = perimeter_all,
             aes(x = per_h,
                 y = per_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = perimeter_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = perimeter_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Herbarium perimeter (mm)",
       y = "Initial perimeter (mm)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 375, by = 50), 
                     limits = c(0, 375)) +
  scale_y_continuous(breaks = seq(from = 0, to = 375, by = 50), 
                     limits = c(0, 375))
perimeter_ih_plot
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 63 rows containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of per_h(x) as Predictor of per_i(y) with a species facet.
perimeter_ih_plot_species <- ggplot() +
  
  geom_point(data = perimeter_all,
             aes(x = per_h,
                 y = per_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = perimeter_ih_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = perimeter_ih_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Herbarium perimeter (mm)",
       y = "Initial perimeter (mm)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  perimeter_ih_plot_species
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

### Surface Area

     #....Scatterplot of surf_h(x) as a Predictor of surf_i(y)....
surf_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = surf_all,
             aes(x = surf_h,
                 y = surf_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = surf_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = surf_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Herbarium Surface Area (mm^2)",
       y = "Initial Surface Area (mm^2)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 450, by = 75), 
                     limits = c(0, 450)) +
  scale_y_continuous(breaks = seq(from = 0, to = 450, by = 75), 
                     limits = c(0, 450))
surf_ih_plot
## Warning: Removed 159 rows containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of surf_h(x) as Predictor of surf_i(y) with a species facet.
surf_ih_plot_species <- ggplot() +
  
  geom_point(data = surf_all,
             aes(x = surf_r,
                 y = surf_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = surf_ih_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = surf_ih_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Herbarium Surface Area (mm^2)",
       y = "Initial Surface Area (mm^2)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  surf_ih_plot_species

### Weight

     #....Scatterplot of weight_h(x) as a Predictor of weight_i(y)....
tdmc_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = tdmc_all,
             aes(x = weight_h,
                 y = weight_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = tdmc_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = tdmc_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Herbarium weight (g)",
       y = "Initial weight (g)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 4850, by = 810), 
                     limits = c(0, 4850)) +
  scale_y_continuous(breaks = seq(from = 0, to = 4850, by = 810), 
                     limits = c(0, 4850))
tdmc_ih_plot
## Warning: Removed 54 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 28430 rows containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of weight_h(x) as Predictor of weight_i(y) with a species facet.
tdmc_ih_plot_species <- ggplot() +
  
  geom_point(data = tdmc_all,
             aes(x = weight_h,
                 y = weight_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = tdmc_ih_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = tdmc_ih_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Herbarium weight (g)",
       y = "Initial weight (g)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


 tdmc_ih_plot_species
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).

### Volume

     #....Scatterplot of vol_h(x) as a Predictor of vol_i(y)....
volume_ih_plot <- ggplot() +
  
  # plot the raw data as points
  geom_point(data = volume_all,
             aes(x = vol_h,
                 y = vol_i),
             alpha = 0.5,
             color = model_col,
             shape = 21) +
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence interval 
  geom_ribbon(data = volume_ih_predictions,
              aes(x = x,
                  y = predicted,
                  ymin = conf.low,
                  ymax = conf.high),
              alpha = 0.2,
              fill = model_col) +
  
  # plot the prediction line
  geom_line(data = volume_ih_predictions,
            aes(x = x,
                y = predicted),
            color = model_col,
            linewidth = 1) +
  
  # labels
  labs(x = "Herbarium volume (mL)",
       y = "Initial volume (mL)") +
  
  # controlling axes to make the plot look square
  # this makes it easier to see the difference between the model and the 1:1 line
  scale_x_continuous(breaks = seq(from = 0, to = 15, by = 2.5), 
                     limits = c(0, 15)) +
  scale_y_continuous(breaks = seq(from = 0, to = 15, by = 2.5), 
                     limits = c(0, 15))
volume_ih_plot
## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 29 rows containing missing values or values outside the scale range
## (`geom_line()`).

  #.Scatterplot of vol_h(x) as Predictor of vol_i(y) with a species facet.
volume_ih_plot_species <- ggplot() +
  
  geom_point(data = volume_all,
             aes(x = vol_h,
                 y = vol_i,
                 color = species),
             alpha = 0.5,
             shape = 21) + 
  
  # plot a 1:1 reference line
  geom_abline(slope = 1, 
              intercept = 0,
              linetype = 2,
              linewidth = 1,
              color = ref_line_col) +
  
  # plot the confidence intervals
  geom_ribbon(data = volume_ih_predictions_species,
            aes(x = x,
                y = predicted,
                ymin = conf.low,
                ymax = conf.high,
                fill = species,
                group = species),
            alpha = 0.2) +
  
  # plot the model predictions
  geom_line(data = volume_ih_predictions_species,
            aes(x = x,
                y = predicted,
                color = species,
                group = species)) +
  
  # labels
  labs(x = "Herbarium volume (mL)",
       y = "Initial volume (mL)") +
  
  # theme
  theme(legend.position = "none",
        strip.background = element_blank(),
        strip.text = element_text(size = 12)) +
  
  # facet
  facet_wrap(~ species, nrow = 2,
             scales = "free")


  volume_ih_plot_species

c. Visualizing summary graphs for rehydrate models

This section includes visualizations that plot multiple trait models on a single graph

Scatterplot faceted by trait

This scatterplot visualizes predictions for height, width, perimeter, weight, surface area, thickness, and volume as a function of inital trait values, faceted by trait.

# surf_ir_predictions, perimeter_ir_predictions, tdmc_ir_predictions, height_ir_predictions, width_ir_predictions, thick_ir_predictions, volume_ir_predictions

     #.........................Data Wrangling.........................

#creating a list of all the predictions for rehydrate traits
init_rehy_df <- rbind(surf_ir_predictions, # rbind stacks data frames on top of each other if they have the same columns
                      perimeter_ir_predictions, 
                      tdmc_ir_predictions, 
                      height_ir_predictions, 
                      width_ir_predictions, 
                      thick_ir_predictions, 
                      volume_ir_predictions) |> 
  drop_na() # not sure why there are NA rows at the end of the data frame, but they come from one of the prediction data frames

     #.........A scatterplot with all seven model predictions.........
ggplot(data = init_rehy_df, # using the data frame created above
       aes(x = x, # x-axis
           y = predicted )) + # y-axis
  geom_abline(slope = 1, # adding a dashed grey 1:1 reference line
              intercept = 0,
              color = "grey",
              linetype = 2) + 
  # geom_ribbon(aes(fill = trait, # adding 95% CI 
  #                 ymin = conf.low,
  #                 ymax = conf.high),
  #             alpha = 0.2) +
  geom_line(aes(color = trait), 
            linewidth = 1) + # plotting model predictions
  labs(x = "Rehydrated value", # renaming the axes
       y = "Initial value") +
  theme(legend.position = "none") +
  facet_wrap(~trait, scales = "free")

Highlighting 1:1 relationships

This scatterplot visualizes predictions for height, width, perimeter, surface area, thickness, and volume as a function of inital trait values. Prediction lines are on the same graph, and models that follow a 1:1 relationship are highlighted.

     #..............................Setup.............................

#Using the dataframe with all our rehydrate trait predictions to assign an opacity
#value to all traits. 3 signifies a close 1:1 relationship while 0.1 is a weak 1:1.
init <- init_rehy_df 
init23 <- init[init$trait != "thallus dry matter content"]
init23$opacity <- ifelse(init23$trait == "height", 3,
                      ifelse(init23$trait == "width", 3,
                             ifelse(init23$trait == "perimeter", 0.1, 
                                    ifelse(init23$trait == "surface area", 0.1, 
                                                  ifelse(init23$trait == "thickness", 0.1, 
                                                         ifelse(init23$trait == "volume", 0.1, NA ))))))
     #........manually offsetting height and width line labels........
label_df <- init23 |> 
  filter(x == max(x)) |> 
  mutate(label = stringr::str_to_title(trait),
         nudge_y = case_when(
           trait == "height" ~ 0.7,
           trait == "width" ~ -0.3,
           TRUE ~ 0
         ))

     #........................generating titles.......................
title <- "Two predictive trait models based on rehydrated values closely follow a 1:1 relationship"
subtitle <- "The traits that follow this pattern are Height and Width"

     #..............assigning each trait a unique color...............
showtext_auto()
trait_palette <- c("surface area" = "#ffbf00", 
                   "perimeter" = "#be8333", 
                  # "thallus dry matter content" = "#54662c", 
                   "height" = "#009bb0", 
                   "width" = "#124c54", 
                   "thickness" = "#368000", 
                   "volume" = "#b4450e")

     #........scatterplot with all rehydrate trait predictions........
     #......rehydrated trait value (x) & initial trait values (y).....
     #..excludes tdmc, and highlights traits with 1:1 relationships...
 plot_ir_1s <- init23 |> 
  ggplot(aes(x, predicted, color = trait, alpha = opacity))+
  geom_line()+
  geom_abline(slope = 1, # adding a dashed grey 1:1 reference line
              intercept = 0,
              color = "grey",
              linetype = 2) + 
  # geom_ribbon(aes(fill = trait, # adding 95% CI (uncomment this if you want to see it)
  #                 ymin = conf.low,
  #                 ymax = conf.high),
  #             alpha = 0.2) +
  geom_line(
            linewidth = 0.5) + # plotting model predictions
  labs(x = "Rehydrated Value", # renaming the axes
       y = "Predicted Initial Value", 
       title = title,
       subtitle = subtitle) +
  scale_color_manual(values = trait_palette)+ #assigning values a color based on trait from our palette
  theme(legend.position = "none", 
        plot.background = element_blank())+
 #formatting line labels
    geom_text(
    data = init23 |> filter(x == max(x)) |> mutate(label = stringr::str_to_title(trait)),
    aes(label = label, alpha = opacity), 
    hjust = 0, 
    nudge_x = 0.1,
    nudge_y = label_df$nudge_y,
    size = 5, 
    family = "lex", 
    #alpha = opacity
    #fontface = "light"
  )+
   #formatting text and margins
  coord_cartesian(xlim = c(0, 7.2), 
                  clip = "off")+
  theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, "cm"))+
   theme(axis.title.x = element_text(family = "lex",
                                 size = 14,
                                 hjust = 0.5,
                                 margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")),
         axis.title.y = element_text(family = "lex",
                                 size = 14,
                                 hjust = 0.5,
                                 margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")), 
          plot.title = element_text(family = "lex", 
                              size = 17, 
                              hjust = 0.5,
                              margin = margin(t = 0, r = 1, b = 0.3, l = 0, "cm")),
    plot.subtitle = element_text(family = "lex",
                                 size = 15,
                                 hjust = 0.5,
                                 margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")))+
  scale_x_continuous(expand = c(0, NA))+
  scale_y_continuous(limits = c(0.8,NA))
 
 plot_ir_1s
## Warning in y + params$y: longer object length is not a multiple of shorter
## object length
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_text()`).

Highlighting relationships that deviate from 1:1

This includes code for a scatterplot of height, width, perimeter, surface area, thickness, and volume on the same graph. Models that deviate from a 1:1 relationship are highlighted.

#.........................assigning opacity values (manually).........................

#traits that follow a 1:1 relationship are assigned 1, while deviants from 1:1
#are assigned 3. 
showtext_auto()
init <- init_rehy_df 
init2 <- as.data.frame(init) |> 
  mutate(trait = as.factor(trait))

init2 <- init[init$trait != "thallus dry matter content"]
init2$opacity <- ifelse(init2$trait == "height", 1,
                      ifelse(init2$trait == "width", 1,
                             ifelse(init2$trait == "perimeter", 3, 
                                    ifelse(init2$trait == "surface area", 3, 
                                                  ifelse(init2$trait == "thickness", 1, 
                                                         ifelse(init2$trait == "volume", 3, NA ))))))

#.........................titles and palettes.........................

title <- "Predictive models based on rehydrated measurements vary by trait"
subtitle <- "The traits that deviate the most are Perimeter, Volume, and Surface Area"

showtext_auto()
trait_palette <- c("surface area" = "#ffbf00", 
                   "perimeter" = "#be8333", 
                  # "thallus dry matter content" = "#54662c", 
                   "height" = "#009bb0", 
                   "width" = "#124c54", 
                   "thickness" = "#368000", 
                   "volume" = "#b4450e")

#.........................plotting.........................
#......rehydrate trait values (x) & initial trait values (y).....

plot_ir_deviants <- init2 |> 
  ggplot(aes(x, predicted, color = trait, alpha = opacity))+
  geom_line()+
  geom_abline(slope = 1, # adding a dashed grey 1:1 reference line
              intercept = 0,
              color = "grey",
              linetype = 2) + 
  # geom_ribbon(aes(fill = trait, # adding 95% CI (uncomment this if you want to see it)
  #                 ymin = conf.low,
  #                 ymax = conf.high),
  #             alpha = 0.2) +
  geom_line(
            linewidth = 1.5) + # plotting model predictions
  labs(x = "Rehydrated Value", # renaming the axes
       y = "Predicted Initial Value", 
       title = title,
       subtitle = subtitle) +
  scale_color_manual(values = trait_palette)+ #assigning palette colors to traits
  theme(legend.position = "none", 
        plot.background = element_blank())+
  geom_text( #formatting model labels
    data = init2 |> filter(x == max(x)) |> mutate(label = stringr::str_to_title(trait)),
    aes(label = label), 
    hjust = 0, 
    nudge_x = 0.1, 
    size = 5, 
    family = "lex"
    #fontface = "light"
  )+
  #manually adjusting text elements and plot dimensions
  coord_cartesian(xlim = c(0, 7.2), 
                  clip = "off")+
  theme(plot.margin = margin(t = 1, r = 1, b = 1, l = 1, "cm"))+
   theme(axis.title.x = element_text(family = "lex",
                                 size = 14,
                                 hjust = 0.5,
                                 margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")),
         axis.title.y = element_text(family = "lex",
                                 size = 14,
                                 hjust = 0.5,
                                 margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")), 
          plot.title = element_text(family = "lex", 
                              size = 17, 
                              hjust = 0.5,
                              margin = margin(t = 0, r = 0, b = 0.3, l = 0, "cm")),
    plot.subtitle = element_text(family = "lex",
                                 size = 15,
                                 hjust = 0.5,
                                 margin = margin(t = 0, r = 0, b = 0.5, l = 0, "cm")))+
  scale_x_continuous(expand = c(0, NA))+
  scale_y_continuous(limits = c(0.8,NA))

plot_ir_deviants
## Warning: Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 6 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_text()`).